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Fluctuation  induced  almost  invariant  sets 


Ira  B.  Schwartz1  and  Lora  Billings2 

1  US  Naval  Research  Laboratory,  Code  6792,  Nonlinear  Systems  Dynamics  Section, 

Plasma  Physics  Division,  Washington,  DC  20375 

2 Montclair  State  University,  Department  of  Mathematical  Sciences,  Upper  Montclair,  NJ  07043 

We  consider  the  approximation  of  fluctuation  induced  almost  invariant  sets  arising  from 
stochastic  dynamical  systems.  We  describe  the  dynamical  evolution  of  densities  via  the  SFP 
operator.  Given  a  stochastic  kernel  with  a  known  distribution,  approximate  almost  invariant 
sets  are  found  by  translating  the  problem  into  an  eigenvalue  problem  derived  from  reversible 
Markov  processes.  Two  examples  of  the  methods  are  used  to  illustrate  the  technique. 

I.  INTRODUCTION 

Transport  in  dynamical  systems  is  an  important  topic  of  active  research  in  that  it  applies  to 
many  problems  in  physical  systems  from  low  to  high  dimensions.  Most  of  the  tools  for  handling 
transport  from  a  geometric  point  view  have  come  from  modern  dynamical  systems,  such  as  [16]. 
The  effects  of  noise  on  nonlinear  dynamical  systems  have  also  become  an  important  topic  in  recent 
years.  Examples  occur  in  noise  induced  instabilities  arising  in  deterministic  stable  dynamics,  such 
as  escape  from  a  potential  [15],  analysis  of  stochastic  bifurcations  including  exploration  of  the 
interaction  of  noise  and  global  bifurcation  based  on  underlying  chaotic  saddles  [11],  noise  induced 
escape  using  the  Hamiltonian  theory  of  large  fluctuations  [6],  the  theory  of  quasi-potentials  [9,  10], 
or  a  variational  formulation  of  optimal  escape  paths  [13].  It  is  well  known  that  noise  can  excite 
unstable  chaotic  structures  while  destroying  regular  periodic  dynamics,  but  most  studies  consider 
noise  induced  chaos  occurring  near  a  bifurcation,  such  as  a  saddle-node  point  or  a  crisis  of  chaotic 
attractors,  which  leaves  a  chaotic  saddle  present. 

As  an  alternative  to  the  above  analysis  methods,  set  theoretic  methods  have  been  used  to 
quantify  the  stochastic  attractors,  as  well  as  transport  in  stochastic  dynamical  systems.  While  for 
continuous  dynamical  systems  with  noise  one  solves  a  master  equation  (such  as  Fokker-Planck)  for 
the  probability  density  function  [15],  discrete  dynamical  systems  benefit  from  using  the  Frobenius- 
Perron  operator  (FP)  formalism  [12]  to  determine  density  evolution.  In  [4],  the  idea  of  examining 
transition  functions  for  deterministic  chaotic  systems  generated  the  machinery  for  discretizing  the 
FP  operator,  and  was  applied  to  kernels  which  correspond  to  deterministic  systems.  Using  similar 
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ideas,  the  FP  operator  was  explicitly  extended  to  stochastic  kernels  in  [1],  where  it  was  applied  to 
probabilistic  transport  in  epidemic  models  and  discrete  low  dimensional  examples.  A  full  treatment 
and  analysis  of  the  discretization  of  the  FP  operator  is  given  in  [3],  as  well  as  transport  from  one 
invariant  region  to  another. 

More  recently,  the  transport  techniques  have  been  extended  to  extract  almost  invariant  sets 
from  dynamical  systems,  sets  which  are  structures  in  which  the  dynamics  may  remain  for  a  long 
period  of  time  prior  to  leaving  the  region.  [8]  The  techniques  make  use  of  results  from  graph  theory 
on  reversible,  or  detailed  balance,  Markov  processes  of  the  probabilistic  transition  dynamics  [7], 
and  we  follow  a  similar  approach  here  for  stochastic  kernels. 

The  layout  of  the  remaining  part  of  the  paper  is  as  follows.  In  Section  2  we  describe  the  general 
theory  for  the  stochastic  FP  operator.  Both  continuous  and  discrete  versions  are  given.  Methods 
for  computing  approximate  almost  invariant  sets  are  described.  Section  3  has  two  examples  of 
fluctuation  induced  almost  invariant  sets.  One  is  a  simple  bi-stable  one  dimensional  map  and  the 
other  is  an  example  from  epidemiology.  Conclusions  are  summarized  in  Section  4. 

II.  GENERAL  THEORY 

A.  Stochastic  Frobenius-Perron  operator 

We  define  the  Frobenius-Perron  (FP)  operator  for  a  deterministic  function  the  following  way. 
Let  F  be  map  for  a  discrete  dynamical  system  acting  on  a  set  M  C  Rn  defined  by 

F  :  M  — >  M,  x  e- >  F{x).  (1) 

We  define  an  associated  dynamical  system  over  the  space  of  densities  of  ensembles  of  initial  condi¬ 
tions, 

Pf:L\M)^L\M),  p(x)  »  PF[p(x)].  (2) 

This  FP  operator  (Pp)  is  defined  by  [12],  [17] 

PF[p(x)\=  [  5(x  —  F(y))p(y)dy,  (3) 

JM 

acting  on  probability  density  functions  p  £  Ll(M). 

Now  consider  the  stochastically  perturbed  dynamical  system  Fv  :  Rn  — >  Rn,  x  e- >  F(x)  +  rj, 
where  r/  is  a  random  variable  having  probability  density  function  i/(x).  The  stochastic  Frobenius- 


3 


Perron  (SFP)  operator  is  defined  to  be 

PF[p(x)\=  u(x  -  F{y))p(y)dy. 

JM 

For  our  applications,  we  will  use  a  normal  distribution  in  Rn  given  by 


v(x)  = 


1 


JItuj- 


exp(— ||.t||2/2(j2) 


(4) 


(5) 


with  mean  ( x )  =  0  and  standard  deviation  a.  For  a  discrete  time  system  with  constantly  applied 
stochastic  perturbation  having  normal  distribution,  it  is  clear  that  the  SFP  operator  is 

l  r 


pf\p(x)\  = 


V2 


2<r2 


7 RTZ  JM 


"  p(y)dy. 


(6) 


B.  The  discrete  theory 


In  order  to  compute  finite  dimensional  approximations  to  the  SFP  operator,  we  approximate 
functions  in  L1(M),  with  discretely  indexed  basis  functions  {^(a:)}?^  C  Ll{M).  The  finite 
dimensional  linear  subspace  is  generated  by  a  subset  of  the  basis  functions,  Ay  =  span({(pl(x)}1^=l  ).i 
such  that  fa  e  L1(M)Vi.  One  approximates  density  p(x)  by  the  finite  sum  of  basis  functions, 
p(x )  —  YliLi  ci(t)i{x)t  where  4>i(x)  =  XbJx),  xb  is  an  indicator  defined  on  boxes  {BJ  covering  M. 
The  infinite-dimensional  matrix  is  approximated  by  the  N  x  N  matrix, 

A(Bi,Bj)  =  Aitj  =  (■ PFA<t>j\Ai )  =  [  PFa[<t>j(x)]Mx)dx  (7) 

JM 

for  1  <  i,j  <  N.  Therefore,  a  transport  matrix  entry  Ay  value  represents  how  mass  flows  from 
cell  Bi  to  cell  Bj.  Normalizing  Ay  by  the  measure  of  Bi  yields  the  probability  that  a  point  in  Bi 
has  its  image  in  Bj ,  which  is  denoted  by  Py .  P  is  therefore  a  stochastic,  or  probability,  transition 
matrix. 

Almost  invariance  can  now  be  considered  in  terms  of  a  given  partition  which  covers  the  set  M. 
We  can  think  of  a  set  being  almost  invariant  if  its  self  transition  A(Pj,  Bi)  is  close  to  1.  Therefore, 
we  may  adopt  the  following  description  of  almost  invariance.  We  say  an  almost  invariant  set 
consists  of  partition  N  =  {B i,  •  •  •  ,6m}  such  that  A(d)  =  4  A(Bi,Bi)  is  maximized  over  all 

partitions.  [8]  In  order  to  compute  the  almost  invariant  sets,  we  will  require  some  modifications  to 
the  transport  matrix. 


C.  Properties  of  the  transport  matrix 

Given  a  stochastic  transition  matrix  P  where  rows  sum  to  unity,  there  are  several  quantities 
that  can  be  easily  computed  from  the  Perron-Frobenius  theorem: 
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1.  Probabilities  evolve  discretely  according  to  7r(n+b  =  P1  RR. 

2.  There  exists  a  left  eigenvector  it  of  P  which  approximates  the  invariant  probability  density 
of  the  dynamical  system.  It  satisfies  7 t1  =  i ttP,  since  it  has  an  eigenvalue  of  1. 

To  consider  almost  invariant  sets  of  Markov  chains,  we  will  need  to  create  a  transition  matrix  that 
is  reversible.  The  idea  of  using  reversible  Markov  processes  to  find  almost  invariant  sets  comes 
from  an  algorithm  in  [7],  which  in  turn  uses  results  derived  from  ergodic  graph  theory  ideas  applied 
to  image  segmentation.  [14] 

Let  P  be  a  stochastic  primitive  matrix.  [18]  For  such  a  matrix,  there  exists  a  unique  probability 
vector,  such  that  it  is  a  left  eigenvector  with  eigenvalue  1;  i.e. ,  n1  =  nTP.  It  also  has  a  right 
eigenvector,  1  =  (1, . . . ,  1)T.  Since  E  ni  =  1>  the  left  and  right  corresponding  eigenvectors  satisfy 
it'1  1  =  1.  The  interpretation  of  7r  is  very  useful,  since  each  component  7Tj  represents  the  fraction 
of  time  the  dynamics  spends  in  box  Rj. 

We  now  consider  the  idea  of  a  reversible  Markov  chain,  or  detailed  balance.  A  reversible  Markov 
chain  is  one  which  satisfies  7 TjPjj  =  TijPji  for  all  i ,  j .  Such  a  condition  implies  7r  is  an  invariant 
distribution  of  P.  A  reversible  Markov  chain  in  equilibrium  is  one  whose  forward  sequence  of 
events  has  the  same  probability  of  the  reverse  sequence,  making  it  difficult  to  state  which  direction 
time  is  going  in  a  real  experiment.  Letting  D  =  diag^^/Hj),  it  is  easy  to  show  detailed  balance  is 
equivalent  to  D2P  =  PT D2 . 

Since  in  general  P  is  not  reversible  as  we  have  defined  it  from  Eq.  7,  we  construct  a  reversible 
Markov  from  our  transition  matrix.  Let  R  =  (P  +  P)/ 2,  where  Rj  =  -^Pji  and  7r  is  an  invariant 
probability  density  of  P.  Then  the  following  properties  may  be  shown  to  be  true: 

1.  D2R  =  RTD2,  so  R  is  reversible. 

2-  E"=i  Rij  =  !• 

3.  7 ttR  =  7 rT. 

4.  If  R  is  reversible,  the  eigenvalues  are  real  and  the  eigenvectors  are  orthogonal. 

Property  4  is  especially  important,  since  it  will  used  to  cluster  the  invariant  sets. 

In  order  to  compute  a  collection  of  sets  which  are  almost  invariant,  one  may  examine  eigenvectors 
of  R,  which  reveal  the  approximate  almost  invariant  partitions.  Our  machinery  based  on  the  SFP 
operator  defined  in  Eq.  4  allows  us  to  do  this  quite  easily.  The  basic  idea  behind  generating  the 
invariant  sets  via  the  eigenstructure  of  R  may  be  understood  as  follows. 
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Invariant  set  A 

Invariant  set 
B 

FIG.  1:  invariant  set  sketch 


Consider  the  deterministic  case  of  Eq.  1  in  the  absence  of  noise,  and  assume  there  exist  two 
basins  of  attraction,  A  and  B.  Each  of  these  sets  are  made  up  of  aggregates  of  the  appropriate 
union  of  sets  from  the  partition  N.  Given  any  non-empty  index  set  /,  define  its  characteristic  vector 
1/  by  1  /  j  =  1  for  all  iel ,  and  zero  otherwise.  Then  the  one  step  transition  probability  from  a  set 
A  to  B  is  given  by 


wn(A,B) 


£ 


aeA,beB 


£ 


aeA 


'Ka.Pab 


<  1b,  PI  A  >7t 
^  1  At  lA  ^ ir 


(8) 


where  <,  >n  is  an  inner  product  weighted  by  the  distribution  7r;  <  x,y  >n=  xTD2y.  (This  should 
be  compared  with  Eq.  7.)  Clearly,  if  A  is  an  invariant  for  the  deterministic  case,  wn(A,A)  =  1. 
Using  the  FP  machinery,  it  is  clear  that  if  partition  N  may  be  decomposed  into  two  disjoint  invariant 
sets,  then  the  stochastic  case  should  be  a  perturbation  of  these  sets  [5].  In  the  deterministic  case, 
if  A  and  B  are  invariant  and  uncoupled,  then  there  exist  right  eigenvectors  of  R  corresponding  to 
each  diagonal  block  having  eigenvalue  of  1,  with  eigenvector  entries  equal  to  ones.  That  is,  there 
exists  an  eigenspace  spanned  by  the  vectors 


XA  =  (lS,0)  and  Xb  =  (0,1b). 


(9) 


Using  Eq.  9,  we  can  write  any  basis  of  the  eigenspace  in  terms  of  eigenvectors  which  are  constant 
on  the  invariant  sets.  Perturbation  theory  states  that  we  expect  this  structure  to  persist  in  the 
stochastic  cases  as  well,  where  now  the  sets  A  and  B  are  almost  invariant  in  the  diagonal  blocks 
in  Fig.  1,  while  the  off  diagonal  blocks  represent  sets  in  which  transition  occur  of  the  form  A  —>  B, 
and  vice  versa.  In  [1],  we  have  exploited  this  idea  using  the  graph  theoretic  notions  presented  in 
detail  in  [3]. 

Once  we  have  the  eigenstructure  of  R,  we  may  associate  different  signs  with  different  eigen- 
structures.  That  is,  we  find  the  first  few  eigenvalues  of  R  clustered  near  unity,  and  examine  the 
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FIG.  2: 

Plots  of  the  second  and  third  eigenvectors  (left  and  right)  of  the  stochastic  one  dimensional  map  given  in 

Eq.  10  with  a  standard  deviation  of  a  =  0.005. 

eigenstates,  which  will  be  a  mixture  of  almost  invariant  sets  formed  from  the  above  basis  elements. 
(In  [7],  they  use  clustering  arguments  based  on  fuzzy  set  theory,  which  we  do  not  do  here.)  Since 
the  first  eigenvector  of  R  is  a  vector  of  ones,  and  Property  4  states  that  the  second  eigenvector  must 
be  orthogonal  to  the  first,  we  expect  the  second  eigenvector  to  see  if  it  contains  almost  invariant 
sets  emanating  from  two  basins  of  attraction.  Therefore,  we  may  expect  to  use  a  signed  structure 
to  associate  the  right  eigenvectors  of  a  reversible  matrix  with  each  almost  invariant  set.  (There 
do  exist  other  approaches,  but  these  will  be  filled  in  the  discussion.)  We  illustrate  this  with  two 
examples,  both  of  which  consist  of  bi-stable  invariant  attractors  in  the  deterministic  cases. 

III.  EXAMPLES 

In  this  section  we  present  two  examples:  one  a  discrete  map,  and  another  a  driven  flow  perturbed 
with  a  discrete  noise  source. 


A.  A  one  dimensional  map 


The  one-dimensional  piecewise  linear  map,  /  :  [0, 1]  — >  [0, 1], 


/(*) 


— 0.9x  +  0.09 
2.5x  -  0.25 

< 

— 2.5x  +  2.25 
-0.9xTl.81 


if  x  <  0.1, 

if  0.1  <  x  <  0.5, 

> 

if  0.5  <  x  <  0.9, 
if  0.9  <  x, 


(10) 


There  are  two  stable  fixed  points,  x  =  9/19  and  x  =  18.1/19,  with  trapping  regions  [0,0.1]  and 
[0.9, 1],  respectively.  There  are  also  two  unstable  fixed  points,  x  =  1/6  and  9/14.  This  is  a  one- 
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FIG.  3: 

Plots  of  the  second  and  third  eigenvectors  (left  and  right)  of  the  stochastic  one  dimensional  map  given  in 
Eq.  10  with  a  standard  deviation  of  a  =  0.05.  Notice  that  the  small  noise  case  in  Fig.  2  picks  up  almost 
invariant  sets  with  more  detail  than  the  larger  noise  case  here. 

dimensional  bi-stable  system.  Note  that  in  this  example,  we  know  that  each  of  the  two  basins  forms 
a  Cantor  set  in  [0, 1],  with  a  naturally  finer  and  finer  scale  associated  with  each  successive  pre-image 
of  the  trapping  region.  Given  such  a  structure,  it  is  clear  that  the  basins  are  intertwined,  as  shown 
in  [3] .  The  ordering  of  the  sets  is  therefore  quite  easy  to  determine  so  that  a  block  diagonal  matrix 
may  be  constructed  for  invariant  measures. 

In  Figs.  2  and  3,  we  have  implemented  the  SFP  operator  in  Eq.  4  with  the  normal  distribution 
on  the  map  given  by  Eq.  10.  Using  the  invariant  distribution,  we  computed  a  reversible  Markov 
process.  The  first  right  eigenvector  was  found  to  be  uniform,  as  expected,  and  the  second  and 
third  eigenvectors  are  plotted.  The  signs  of  the  entries  in  the  second  eigenvectors  denote  the 
almost  invariant  sets  for  two  distinct  noise  sources.  Figure  2  shows  the  eigenvectors  corresponding 
to  a  standard  deviation  of  a  =  0.005,  and  Fig.  3  has  a  higher  standard  deviation  of  a  =  0.05. 
One  interesting  observation  in  the  second  eigenvector  for  smaller  noise  is  that  the  sets  which  are 
the  union  of  intervals  which  are  positive  appear  to  come  from  the  attractor  on  the  left,  while  the 
negative  values  come  from  the  attractor  on  the  right.  In  the  large  noise  case,  it  appears  that  the 
almost  invariant  set  becomes  smoother,  since  the  standard  deviation  used  in  the  SFP  kernel  is  an 
order  of  magnitude  larger. 

B.  Predicting  almost  invariant  sets  from  epidemics 

We  consider  a  modified  epidemic  model  for  the  spread  of  disease  where  the  state  variables  are 
susceptibles  S  and  infectives  I.  (See  [2]  for  a  derivation  and  notation  of  parameters.)  The  modified 


SI  model  (MSI),  given  by: 


S'{t)  =  n~  n S(t )  -  /3{t)I{t)S(t) 

m  =  (^)  mms(t)  -^+ a)m.  (n) 

(3(t)  =  (3 o(l  +  5cos2irt ) 

The  parameter  values  are  set  to  the  standard  used  for  measles,  and  are  fixed  at  //  =  0.02,  a  = 
1/0.0279,  7  =  1/0.01,  p0  =  1575.0,  6  =  0.1. 

The  solutions  of  the  MSI  model  as  well  as  the  bifurcation  diagrams  for  the  above  set  of  param¬ 
eters  agree  quite  well  with  those  of  the  full  SEIR  model.  For  a  full  description  of  the  deterministic 
solutions,  see  [2].  Since  the  MSI  model  is  periodically  driven  with  period  one  and  both  S  and  I 
are  fractions  of  the  population,  it  may  be  a  viewed  as  a  two-dimensional  map  of  the  unit  box  into 
itself.  The  stochastic  model  is  considered  to  be  discrete  as  well  for  the  purposes  of  this  paper. 
That  is,  noise  is  added  to  the  population  rate  equations  periodically  (period=l)  at  the  same  phase 
having  mean  zero  and  standard  deviation  a.  The  dynamics  may  then  be  represented  as  a  map, 
F  :  3ft2  — ►  SR2,  where 

(S,I)  (t  +  l)  =  F[(S,I)  (t)]+V(t). 

Here  77(f)  is  a  discrete  noise  term. 

In  [2] ,  it  was  shown  that  the  existence  of  bi- instability  is  a  sufficient  condition  for  noise  to  excite 
chaos.  That  is,  there  exists  a  sufficiently  large  standard  deviation  in  which  new  unstable  orbits 
are  created  and  sampled.  For  the  parameters  listed,  in  the  zero  noise  case  there  exist  period  one 
and  period  three  attractors.  The  transport  due  to  the  addition  of  noise  using  Eq.  4  was  originally 
done  in  [1,  3].  Here  we  extend  the  results  to  extract  the  almost  invariant  sets  by  analyzing  the 
eigenstructure  of  the  second  and  third  vectors  of  R. 

We  show  the  small  noise  almost  invariant  sets  (a  =  0.005)  in  Fig.  4,  while  larger  noise  (a  = 
0.035)  results  are  seen  in  Fig.  5.  Notice  that  in  the  left  panel  of  Fig.  4,  there  are  two  almost  invariant 
sets  with  opposite  signs.  The  red  regions  correspond  to  sets  which  are  near  the  deterministic  basin 
boundary  for  the  period  3  attractor,  while  the  dark  blue  correspond  to  sets  near  a  period  2  attractor. 
The  third  eigenvector  begins  to  show  something  new.  The  red  and  yellow  regions  correspond  to 
the  almost  invariant  sets  contained  in  the  second  eigenvector,  but  the  dark  blue  region  points  to  an 
almost  invariant  set  of  an  unstable  region.  This  region  is  what  we  would  call  a  fluctuation  induced 
almost  invariant  set. 

Further  examination  of  the  fluctuation  induced  almost  invariant  set  may  be  made  by  considering 
a  larger  noise  amplitude.  Intuitively,  we  expect  the  noise  to  allow  the  dynamics  to  sample  larger 
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FIG.  4:  The  plots  show  a  projection  of  the  second  and  third  eigenvectors  of  the  reversible  matrix  R  for  the 
MSI  model  in  Eq.  11.  The  parameters  are  listed  in  the  text.  The  standard  deviation  is  a  =  0.005,  which  is 
a  small  noise  case. 
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FIG.  5:  Same  as  in  Fig.  4,  but  with  large  noise  standard  deviation,  a  =  0.035. 


regions  of  phase  space.  The  current  machinery  allows  us  to  describe  more  fully  how  these  regions 
behave.  Figure  5  clearly  shows  that  the  effects  of  noise  enlarge  the  almost  invariant  regions.  In  the 
left  panel,  the  period  3  regions  are  still  distinct,  but  now  the  period  two  region  looks  like  a  large 
region  centered  around  a  period  unstable  point.  In  the  right  panel,  it  is  clear  that  the  fluctuation 
induced  region  is  enlarged  considerably,  and  it  is  expected  that  the  dynamics  spends  a  good  deal 
of  time  in  this  region. 


IV.  CONCLUSIONS 


We  considered  the  approximation  of  fluctuation  induced  almost  invariant  sets  arising  from 
stochastic  dynamical  systems.  We  showed  that  by  starting  with  the  stochastic  FP  operator,  a 
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discrete  transition  matrix  may  be  computed  in  one  step,  the  problem  of  formulating  an  approxima¬ 
tion  to  the  almost  invariant  sets  could  be  translated  into  finding  the  eigenstructure  of  a  reversible 
Markov  process.  The  actual  methods  for  locating  clusters,  or  aggregates  of  sets,  have  its  roots 
in  many  problems,  and  therefore,  there  are  many  alternative  techniques.  The  methods  used  in 
[7,  8]  find  their  origins  in  image  segmentation.  [14]  In  those  works,  as  here,  we  based  our  cluster 
analysis  on  signed  eigenstates  of  a  reversible  Markov  process.  However,  other  techniques  such  as 
fuzzy  set  theory,  greedy  algorithms,  and  graph  theory  are  all  methods  which  may  be  implemented 
to  formulate  a  cluster  analysis  algorithm. 
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